#!/usr/bin/env python3
# coding: utf-8
"""
This module provides a subclass of the pymatgen.core.structure.Molecule class
in order to add new attributes dedicated to the forcefield factory.
"""
import numpy as np
import pymatgen
__all__ = ["Molecule"]
[docs]class Molecule(pymatgen.core.structure.Molecule):
"""
A molecule object inherited from the [pymatgen package](http://pymatgen.org)
with new attributes dedicated to the force field factory.
"""
def __init__(self, species, coords, charge=0,
spin_multiplicity=None, validate_proximity=False,
site_properties=None, bond_orders=None, at_charges=None,
hessian=None, grad=None):
"""
Creates a MutableMolecule.
Args:
species: list of atomic species. Possible kinds of input include a
list of dict of elements/species and occupancies, a List of
elements/specie specified as actual Element/Specie, Strings
("Fe", "Fe2+") or atomic numbers (1,56).
coords (3x1 array): list of cartesian coordinates of each species.
charge (float): Charge for the molecule. Defaults to 0.
spin_multiplicity (int): Spin multiplicity for molecule.
Defaults to None, which means that the spin multiplicity is
set to 1 if the molecule has no unpaired electrons and to 2
if there are unpaired electrons.
validate_proximity (bool): Whether to check if there are sites
that are less than 1 Ang apart. Defaults to False.
site_properties (dict): Properties associated with the sites as
a dict of sequences, e.g., {"magmom":[5,5,5,5]}. The
sequences have to be the same length as the atomic species
and fractional_coords. Defaults to None for no properties.
at_charges (list): atomic charges of sites in the molecule
bond_orders (dict): Dict of atom pairs with a bond order value.
hessian (3N x 3N array): The hessian matrix of the molecule
grad (3N x 3 array): The (last) gradient of the molecule
"""
super(Molecule, self).__init__(species,
coords, charge=charge,
spin_multiplicity=spin_multiplicity,
validate_proximity=validate_proximity,
site_properties=site_properties)
self.at_charges = at_charges
self.bond_orders = bond_orders
self._hessian = None
self.hessian = hessian
self.grad = grad
[docs] def get_bo_value(self, pair):
""" return the bond order value from the atom index """
iat, jat = pair
if (iat, jat) in self.bond_orders:
return self.bond_orders[(iat, jat)]
elif (jat, iat) in self.bond_orders:
return self.bond_orders[(iat, jat)]
else:
raise KeyError("The bonds (%d, %d) or (%d, %d) do not exist" %
(iat, jat, jat, iat))
[docs] def get_atoms_in_bond(self, pair):
""" return the bond order value from the atom index """
iat, jat = pair
if (iat, jat) in self.bond_orders or (jat, iat) in self.bond_orders:
return self[iat], self[jat]
else:
raise KeyError("The bonds (%d, %d) or (%d, %d) do not exist" %
(iat, jat, jat, iat))
def _get_hessian(self):
return self._hessian
def _set_hessian(self, matrix):
""" set the hessian if the array shape is consistant with the molecule """
nat = len(self)
hessian = np.array(matrix, dtype=np.float)
if (3 * nat, 3 * nat) == hessian.shape:
self._hessian = hessian
else:
raise ValueError("Hessian shape is " + str(hessian.shape)
+ " while the molecule contains %d atoms." % nat
+ "The shape has to be (%d, %d)" % (3 * nat, 3 * nat))
hessian = property(_get_hessian, _set_hessian)